
Info << "Reading field beta" << nl << endl;
volScalarField beta
(
    IOobject
    (
        "beta",
        runTime.timeName(),
        mesh,
        IOobject::MUST_READ,
        IOobject::AUTO_WRITE
    ),
    mesh
);

Info << "Reading field U" << nl << endl;
volVectorField U
(
    IOobject
    (
        "U",
        runTime.timeName(),
        mesh,
        IOobject::MUST_READ,
        IOobject::AUTO_WRITE
    ),
    mesh
);

forAll(beta,cellI)
{
    beta[cellI] = 0.0;
    doubleScalar xx=mesh.C()[cellI].x();
    doubleScalar yy=mesh.C()[cellI].y();

    doubleScalar r0 = 0.15;

    doubleScalar disk_r = std::sqrt(std::pow((xx - 0.5), 2) + std::pow((yy - 0.75), 2));
    doubleScalar cone_r = std::sqrt(std::pow((xx - 0.5), 2) + std::pow((yy - 0.25), 2));
    doubleScalar hump_r = std::sqrt(std::pow((xx - 0.25), 2) + std::pow((yy - 0.5), 2));

    if(disk_r<=r0)
    {
        if(std::abs(xx-0.5)>0.025 || yy>0.85)
        {
            beta[cellI] += 1.0;
        }
    }

    if(cone_r<=r0)
    {
        beta[cellI] += 1 - cone_r / r0;
    }

    if(hump_r<=r0)
    {
        beta[cellI] += (1 + std::cos(M_PI * std::min(hump_r / r0, 1.0))) / 4.0;
    }


    U[cellI] = vector(0.5 - yy, xx - 0.5, 0.0);
}
volScalarField beta_exact(beta);

#include "createPhi.H"

beta.correctBoundaryConditions();
beta.write();
U.write();

doubleScalar L1 = 0.0, Lu1 = 0.0, Ll1 = 0.0;
doubleScalar L2 = 0.0, Lu2 = 0.0, Ll2 = 0.0;

   
